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ABSTRACT 


Our Galaxy and the nearby Andromeda galaxy (M31) are the most massive members of the Local Group, and they seem to be a bound 
pair, despite the uncertainties on the relative motion of the two galaxies. A number of studies have shown that the two galaxies will 
likely undergo a close approach in the next 4—5 Gyr. We used direct N-body simulations to model this interaction to shed light on the 
future of the Milky Way - Andromeda system and for the first time explore the fate of the two supermassive black holes (SMBHs) that 
are located at their centers. We investigated how the uncertainties on the relative motion of the two galaxies, linked with the initial 
velocities and the density of the diffuse environment in which they move, affect the estimate of the time they need to merge and form 
“Milkomeda”. After the galaxy merger, we follow the evolution of their two SMBHs up to their close pairing and fusion. Upon the 
fiducial set of parameters, we find that Milky Way and Andromeda will have their closest approach in the next 4.3 Gyr and merge over 
a span of 10 Gyr. Although the time of the first encounter is consistent with other predictions, we find that the merger occurs later than 
previously estimated. We also show that the two SMBHs will spiral in the inner region of Milkomeda and coalesce in less than 16.6 
Myr after the merger of the two galaxies. Finally, we evaluate the gravitational-wave emission caused by the inspiral of the SMBHs, 
and we discuss the detectability of similar SMBH mergers in the nearby Universe (z < 2) through next-generation gravitational-wave 
detectors. 
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1. Introduction 


The Milky Way (MW) and the Andromeda galaxy (M31) are the 
two main members of the Local Group, which contains more 
than 80 galaxies and has a total mass of roughly 3 — 5 x 10!? Mo 
— (Gonzalez et al. 2014; van der Marel et al. 2012a). The future 
° evolution of the Local Group is essentially driven by the dynam- 
ics of our Galaxy and M31, and it can be considered a promis- 
ing study object to investigate the processes of galaxy formation 
and evolution. Although the physical and dynamical properties 
e e of the MW-M31 system are rather uncertain, it is likely that the 
: > Local Group is gravitationally bound and decoupled from the 
>< general cosmic expansion, and also that the two galaxies will 
«_) not escape the collision and the final merger. However, the time 
C3 at which this merger will occur is still a matter of debate. The 
main purpose of this work and of our previous studies (Schiavi 
et al. 2019a,b) is to shed light on this topic. 
According to some previous simulations (Dubinski et al. 1996; 
Cox & Loeb 2008; van der Marel et al. 2012b), the first close ap- 
proach will likely occur in < 4 Gyr, even though the two galax- 
ies have different initial conditions in all the cited works. Using 
a more recent estimation of the proper motion of M31, van der 
Marel et al. (2019) have obtained a time for the first approach 
equal to 4.5 Gyr. Almost the same result can also be obtained 
when the so-called “timing argument” is employed, which was 
introduced in the pioneering work by Kahn & Woltjer (1959). 
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In the timing argument, MW and M31 are considered as point 
masses on a radial orbit: they started their motion at the Big 
Bang, and after decoupling from the Hubble flow, began to ap- 
proach one another. Even though the timing argument allowed 
obtaining an estimate of the total mass of the Local Group that is 
compatible with the estimate obtained with other methods (e.g., 
Klypin et al. 2002; Widrow & Dubinski 2005), it is unable to 
take the complexity of the dynamics of the galaxy interaction 
into account. The time needed for the completion of the whole 
merger process is highly sensitive not only to the masses of the 
two galaxies, but also to their proper motion and to the density 
of the intergalactic medium (IGM) in which they move. All the 
estimates of the mass of the two galaxies are affected by a rather 
high level of uncertainty. This is due mainly to the presence of 
extended dark matter halos. We have chosen to adopt the values 
estimated in Klypin et al. (2002), defined as the virial masses at 
radius 7299, where the galactic density is 200 times the critical 
density po ~ 1.0 x 10-*° kg/m? (according to the measurements 
of the Hubble constant Hp by Huang et al. (2020) and Aghanim 
et al. (2018)): Myw = 1.0 10! Mo and My3; = 1.6 x 10!? Mo. 
Another source of uncertainty is our poor knowledge of the ac- 
tual size of the two galactic halos. The radial extent of the halo in 
equilibrium models of galaxies developed by Kuijken & Dubin- 
ski (1995) is in the range of 21 — 73 disk scale lengths. The ratio 
of the dark halo virial radius and the galaxy effective radius fall 
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in the same range in the studies by Jiang et al. (2019), Somerville 
et al. (2017), Huang et al. (2017), and Kravtsov (2013). How- 
ever, there is evidence that some of the gaseous circumgalactic 
medium (CGM) extends to even larger distances (Shull 2014). In 
other words, we do not know with sufficient precision where a 
galaxy actually ends, and in our specific case, whether the Milky 
Way and Andromeda are already partially overlapping or if they 
are still well separated. For this reason, as discussed in Section 
2, we have decided to set the halo cutoff radii of the two galaxies 
at the respective tidal radii. Moreover, in Section 5.1 we demon- 
strate that the time of the merger does not depend on galactic 
halos that are more extended than 80 disk scale lengths. Be- 
cause it is evident that the edge of each galaxy gradually fades in 
the IGM, we cannot ignore the effect of this diffuse medium in 
studying the interaction. The density of the IGM is known to be 
four to six times the critical density pọ (Chamaraux & Tadokoro 
1971; Cox & Loeb 2008), but by performing several simulations, 
we obtained that even a small variation in this parameter could 
affect the merger time substantially. 

Our knowledge of the proper motion of M31 relative to us is 
mainly obtained through redshift measurement. This gained us 
an accurate estimate of the only radial component of the rela- 
tive velocity vector of M31: V, ~ 120 km/s (Binney & Tremaine 
1987). The tangential component has been inferred by studying 
the motion of the satellite galaxies of M31 (Loeb et al. 2005; 
van der Marel & Guhathakurta 2008) or by the Hubble Space 
Telescope (HST) and GAIA observations of sources behind M31 
(van der Marel et al. 2012b, 2014, 2016). These estimates span 
from a minimum of V, ~ 17 km/s (van der Marel et al. 2012a) to 
a maximum of V, x 164 km/s (Salomon et al. 2016). The most 
recent estimate is that by van der Marel et al. (2019), who have 
used GAIA DR2 to obtain V, = ST? km/s. We referred to this 
ultimate measurement to fix the orientation of the relative ve- 
locity vector and its radial component (V, = —115.7 km/s). We 
discuss the relation between the initial tangential velocity and 
the time of the merger in Section 5.1. 

During the interaction at large scales, we are interested in fol- 
lowing the motion of the two SMBHs in the centers of the two 
galaxies. It is well known that a compact object of mass Myw = 
4.31 x 10° Mo (Gillessen et al. 2009), called SgrA*, is placed at 
the center of the Milky Way. Even though the nucleus of M31 
seems to have a double or triple structure, there is a high proba- 
bility that it might host an SMBH of mass Mi; = 1.4 x 108 Mo 
(Bender et al. 2005). After the merger of the two host galaxies, 
their SMBHs are expected to form a binary that will shrink over 
time through gravitational encounters with field stars. We first 
explore the future evolution of some of the nearest SMBHs and 
their eventual coalescence in the nucleus of the galactic merger 
remnant. In Section 5.2 we discuss the time required for the 
SMBHs to merge and the amount of energy radiated through 
gravitational waves (GWs). 

One of the most effective ways to model galaxy interactions is 
the integration of the N-body problem. While tree codes can sim- 
ulate a large number of collisionless particles in galaxy merger 
simulations, a collisional direct summation N-body code with 
fewer particles is required in this study to follow the SMBH dy- 
namics. Direct N-body codes are highly reliable but computa- 
tionally expensive: this clearly places a limit on the number of 
particles involved in the simulations, and prevents us from re- 
solving large and small scales at the same time with good ac- 
curacy. For this reason, as we discuss in Section 3, we chose to 
split the whole study into two parts: in the first, we simulate the 
galaxy interaction at large scales, and in the second, we focus on 
the analysis of the orbital decay of the SMBH binary. 
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2. Galactic model and initial conditions 


Our galaxies were modeled by combining three different com- 
ponents: an exponential disk, a spherical bulge, and a halo, the 
latter two with a Hernquist (Hernquist 1990) density profile. We 
combined these three components with the command magalie 
in the NEMO code (Teuben 1995), which guarantees the stability 
of the whole system. 

The two galaxies have the structure presented in Klypin et al. 
(2002) and in Widrow & Dubinski (2005), that is, nearly the 
same as was used by Cox & Loeb (2008). The main structure 
parameters are summarized in Table 1. 


Milky Way Andromeda 


Scale radius of the disk (kpc) 3.5 5.7 
Core radius of the bulge 0.2 0.2 
Core radius of the halo 3.0 3.0 
Cutoff radius of the halo 98.3 76.5 
Mass of the disk (Mo) 40x100 70x100 
Mass of the bulge 0.2 0.3 
Mass of the halo 23.8 21.6 
Total mass (Mo) 1.0 x 10” 1.6 x 10” 


Table 1. Values of characteristic parameters used in our simulations. 
Where it is not specified, the lengths are in units of the scale radius of 
the disk R4 and masses are in units of the mass of the disk M4. 


The cutoff radii of the halos were chosen as the tidal radii of 
the two galaxies, computed in the point-mass approximation: 


Ft MW 


Muw -( 
Mya 


2 
and rm31 = R — f, mw, d) 
R-rimw 


where r, mw and r, m31 are the two tidal radii and R is the current 
separation between the two galaxies. A snapshot of our galactic 
model is shown in Fig. 1. 

We placed an SMBH in the center of each galaxy, which in the 
first part of our study was modeled as a particle with a mass of 
0.001 times the mass of the whole galaxy. This ratio implies that 
the mass of our SMBHs is significantly higher than the observed 
masses, but this is not relevant for the dynamics of the two galax- 
ies until merging because the two SMBHs are essentially passive 
guests of the hosting galaxies at this phase. We therefore made 
this choice because it was the best setting allowed by our numer- 
ical resolution. We used a number of particles not greater than 
N = 2.6 x 10° , and this constrains the mass of the single parti- 
cle. An SMBH mass of one thousandth of the mass of the galaxy 
is therefore a good compromise between the properties of the 
galaxies and the comparison with an ordinary particle. However, 
during the simulation of the collision at large scales, the two par- 
ticles that represent the SMBHs only have the purpose of better 
identifying the two galactic centers and of observing their rela- 
tive distance at the end of the merger process. 

The Milky Way and Andromeda start to interact at the cur- 
rent distance of 780 kpc (McConnachie et al. 2005; Ribas 
et al. 2005), and their spin vectors are oriented at (0°; —90°) 
and (240°; —30°) in Galactic coordinates, respectively (Dubin- 
ski et al. 1996; Raychaudury & Lynden-Bell 1989). To better 
display the dynamics of the galaxy binary system, we chose a 
reference frame where the x-y plane coincides with the plane of 
the motion. The initial configuration of the two galaxies is shown 
in Fig. 2. 
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Fig. 1. Our model of the Milky Way. The three components (disk, bulge, 
and halo) are shown in different colors. The lower panel is a zoom into 
the innermost region. 
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Fig. 2. Diagram of the initial configuration of the two galaxies in the 
chosen reference frame. 


3. Methods 


Owing to the computational complexity of the simulations, we 
decided to divide the study into two parts. The first examines the 
galaxy interaction at large scales and has the main purpose of 
determining the true time required to form Milkomeda and its 
final density profile. In this section we use a number of particles 
of N = 2.6 x 10° for the simulations with the fiducial set of pa- 
rameters and N = 6.5 x 10* for all the others. The numerical 
integration of the N-body equations of motion was implemented 
with the HiGPUs code (Capuzzo Dolcetta et al. 2013). This pro- 


gram is based on a sixth-order Hermite integration scheme with 
block time steps and directly computes the mutual force between 
each pair of particles, exploiting a parallelization of CPUs and 
GPUs. Owing to the high performance of the HiGPUs code, we 
repeated the simulation several times, changing two parameters 
that were linked with the initial conditions and the external envi- 
ronment: the tangential component of the initial relative velocity 
V, , and the density of the IGM p. This allowed us to investigate 
the correlation between the time of the merger and the galactic 
dynamical properties, together with the effect of the dynamical 
friction exerted by the surrounding diffuse medium. 

In the second part we further study the evolution of the SMBH 
binary that formed after the galaxy merger. Taking the simula- 
tion with the highest resolution and the fiducial set of parame- 
ters, we obtained the Milkomeda density profile and the velocity 
dispersion and modeled the galactic center as an analytic distri- 
bution of matter around the SMBH binary. To simulate the evolu- 
tion of the binary, we used a modified version of the ARWV code 
(Mikkola & Tanikawa 1999; Mikkola & Merritt 2008), which 
integrates the equations of motion taking in account the effect 
of the dynamical friction exerted by a diffuse background dur- 
ing the first phases of the orbital evolution, the post-Newtonian 
(PN) terms when the binary shrinks enough to reach the GW 
emission regime, and the effect of the spins of the merging ob- 
jects (Arca-Sedda & Capuzzo-Dolcetta 2019; Chassonnery et al. 
2019). To connect the new simulation to the previous one, the 
SMBHs orbit was reproduced with the same geometry as found 
at the galaxy merger, while the masses of the two objects were 
set to those known from observational evidences. Through this 
method, we can study the orbital decay of our binary down to 
small spatial scales and infer the coalescence time, the evolution 
of the semimajor axis, the eccentricity, and the power emitted in 
the form of GW. 


4. Intergalactic medium and dynamical friction 


The presence of the IGM affects the time of the galaxy interac- 
tion through the extraction of orbital energy and angular momen- 
tum. We used HiGPUs to simulate the dynamics of the galaxy 
collision in different environments. To take the effect of the IGM 
into account, we modified the HiGPUs code by adding a dy- 
namical friction term in the equation of motion of each particle 
according to the Chandrasekhar formula (Binney & Tremaine 
1987), 


dri N Gm (r; = r;) 
R` 3/2 
a j+i (e + \rj = rf) 
4nG?pM ln A 2X _y 
-PETE fero- ely, 2) 
vè yr 


C 


with X = V,/ V2ø, where o is the IGM velocity dispersion. 
As usual for N-body codes, we introduced the softening param- 
eter € to avoid the divergence of the Newton term at small dis- 
tances: it was fixed at € = 500 pc for ordinary particles and 
e = 50 pc for the two black holes. In the Chandrasekhar term 
we considered p as the density of the IGM, and M and V, as the 
mass and velocity of the galactic core that the particle belongs to. 
Unlike the classical Chandrasekhar formula, which describes the 
effect of the dynamical friction on each star, depending on the 
stellar mass and velocity, in our case, each particle, which has 
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the same mass m, feels the same frictional force as all the oth- 
ers. This force changes with time during the galaxy interaction, 
but at any moment, is the same for every particle. This choice is 
suggested by the need of describing the collective effect of the 
friction on the motion of each galaxy as a whole. 

In all our simulations we fixed the Coulomb logarithm at In A = 
5 and the velocity dispersion of the medium at o = 86.2 km/s, 
obtained from the equipartition of energy for a diffuse medium 
at a temperature of T = 3 x 10° K (Cox & Loeb 2008). 

We compared the case with no IGM with three cases with differ- 
ent values of p: 1.0x10~7° kg/m?, the same as the critical density, 
4.0x 107% kg/m, which is the value estimated by Chamaraux & 
Tadokoro (1971), and 1.0 x 10~?° kg/ m?, about 10 times the crit- 
ical density. As we show in Section 5.1, the time of the merger 
significantly changes for different IGM densities, especially for 
high initial velocities. 


5. Results 
5.1. Galaxy merger 


For all initial conditions, the merger remnant Milkomeda re- 
sembles a giant elliptical galaxy with a density profile similar 
to those of the original two galaxies, as is shown in Fig.3 for 
V, = 57 km/s and p = 4.0 x 107” kg/m?. 


ao 
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Milkomeda 


0.01 0.1 1 10 100 


r (kpc) 


Fig. 3. Density profile of Milkomeda at time t = 10.2 Gyr in the case 
with V, = 57 km/s and p = 4.0 x 107% kg/m’, and those of the Milky 
Way and M31 at t = 0. 


We obtain the time of the merger from the time evolution 
of the distance between the centers of mass. The time of the 
merger is defined here as the time at which the separation is 
0.5% of its initial value. 

Before fixing the outermost edge of our galaxies at the respec- 
tive tidal radii, we investigated the dependence of the time of 
the merger Tm on the cutoff radius R, of the galactic halos. In 
the top panel of Fig. 4 we show that as expected, Tm is strongly 
dependent on the galaxy extension only for low values of Rp. 
For R, > 80Rz , the timing of the process is no longer sensitive 
to this parameter: from this value on, the two halos cover the 
entire distance between the two galaxies. 

We also found that when V, increases, Tm rapidly increases, as 
is shown in the lower panel of Fig.4 in the case of no IGM. 
It is interesting to note that there seems to be a quite accurate 
relation between T,,, and V,: our best fit is Tn œ V,’. 
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Fig. 4. Top: Time of the merger as a function of the cutoff radius of the 
halo R;,, given in units of the scale radius of the disk R4 and assumed 
equal for the two galaxies. The initial tangential velocity here is fixed 
at 50 km/s. Bottom: Correlation between the time of the merger and the 
initial tangential velocity in the case of no IGM. The cutoff radii of the 
halos here are fixed at 70 Ry. 


Repeating the simulation for different values of the IGM den- 
sity, we obtained that the presence of a diffuse medium speeds 
the galaxy interaction up, especially in the case of large V;. In 
Fig. 5 we plot the evolution of the distance between the two cen- 
ters of mass for four different values of p in the case of V, = 57 
km/s (top panel) and the dependence of the time of the merger on 
the IGM density for three different values of V, (bottom panel). 
Even though the time of the merger can significantly change 
when V, or p are varied, we note that the time of the first ap- 
proach is almost constrained in the interval 4—5 Gyr. This means 
that the first part of the galaxy motion is nearly Keplerian be- 
cause the orbital energy dissipation due to the friction exerted by 
the IGM is still not very efficient. The time of the first approach 
is very close to that obtained in the case of a pure radial fall of 
two point masses starting at a distance of 780 kpc with a relative 
radial velocity of -115.7 km/s. After the first encounter, the IGM 
density instead plays a relevant role in the time for the comple- 
tion of the merger. This is mainly due to the enhanced speed of 
the two galaxies at the pericenter. 

Among the ensemble of the performed simulations, we refer to 
the simulation with V, = 57 km/s and p = 4.0 x 107%” kg/m? 
as a fiducial case. According to our analysis, the Milky Way and 
Andromeda will reach their first approach in 4.3 Gyr and will 
merge in 10.0 Gyr. The time for the first pericenter is close to 
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the ~ 4.5 Gyr found by van der Marel et al. (2019), but they did 
not report any value for the time of the final merger. Cox & Loeb 
(2008) obtained the first approach at ~ 2.8 Gyr and the merger 
at ~ 5.4 Gyr. However, we have to consider that they started the 
simulations of the MW-M31 interaction 5 Gyr in the past and 
reached a current transverse velocity that is very different to the 
recently measured velocity. Moreover, they used an IGM density 
that is slightly greater than we considered in our fiducial model. 
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Fig. 5. Top: Separation between the two galaxies as a function of time 
for three different values of IGM density, in the case of V, = 57 km/s. 
Bottom: Time of the merger as a function of the IGM density for three 
values of V,. 


5.2. SMBH merger 


We found that the distance between the two SMBHs located at 
the galactic centers evolves in time, as was previously shown 
for the two centers of mass. The only difference is that after the 
galaxy merger, the SMBH binary stalls at the same fixed dis- 
tance, independently of their initial velocity, as is shown in the 
Fig. 6. This confirms the idea that the orbital decay of the bi- 
nary in the first phase simply follows the dynamics of the two 
stellar systems, but it later depends on the gravitational encoun- 
ters between the binary and the stars orbiting close to the galac- 
tic center. When the volume around the binary is depleted and 
the binary orbit contains a total mass equal to or lower than the 
combined mass of the two SMBHs, the binary stalls. Therefore, 
as expected, this interesting behavior is a function of the num- 
ber of particles that is involved in the simulation: the greater 


the number of particles, the smaller the stalling distance. Nev- 
ertheless, we note that as the numerical resolution increases at 
N > 5.0 x 10%, the stalling distance reaches an asymptotic value 
of ~ 100 pc, that is, about twice the softening parameter of the 
two SMBHs. Fig.7 shows this behavior: for three different val- 
ues of N > 5.0 x 104, the stalling distance does not change. 
This might be the signal that we have reached the lower limit of 
the stalling distance that is allowed by our computational power. 
The density of stars around the binary rapidly drops to zero be- 
cause of the low numerical resolution, and this makes the en- 
ergy loss by dynamical friction inefficient. However, because the 
main purpose of this first simulation is to reproduce the galaxy 
interaction, we cannot expect to simultaneously correctly resolve 
the dynamics at small scales. The stall shown in Fig. 7 is there- 
fore a direct effect of the low spatial resolution and an indirect 
effect of the sampling. 
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Fig. 6. Distance between the two SMBHs as a function of time for three 
different initial tangential velocities. 
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Fig. 7. Distance between the two SMBHs as a function of time for three 
different numerical resolutions of the simulation. 


5.2.1. Orbital decay of the SMBH binary 


To follow the orbital decay of the SMBH binary, we used an- 
other code that improves the treatment of dynamical friction for 
the second part of our study. We chose the simulation with the 
largest number of particles (N = 2.6 x 10°) and the fiducial set 
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of parameters (V, = 57 km/s, p = 4.0 x 10-76 kg/m?) and cal- 
culated the density profile of Milkomeda at time t = 10.2 Gyr, 
soon after the merger, when the system has stabilized and the 
SMBH binary has formed. This last point occurs when the two 
objects approximately reach the so-called influence radius, de- 
fined in Merritt et al. (2007) as the radius of a sphere around the 
two SMBHs that contains twice the sum of their masses. In our 
case, this corresponds to r} % 156.5 pc. 

We used a Dehnen power law (Dehnen 1993), with y = 0.8, 
scale length 50 kpc, and total mass 9.3 x 10! Mo, to fit the 
innermost region of the density profile of Milkomeda and re- 
produce the galactic center as an analytic external potential in 
the ARWV code. The velocity dispersion in the innermost 500 
parsecs (0 = 203.5 km/s) was obtained from the outputs of 
HiGPUs code. 

The values of the masses of the two SMBHs were set to the 
proper values Myw = 4.31 x 10° Mo, and Mi = 14x 108 Mo 
(Gillessen et al. 2009; Bender et al. 2005), keeping as initial con- 
ditions those coming from the last computed orbit of the SMBH 
pair. The orbital integration with the ARWV code was performed 
in a reference frame in which the x-y plane coincides with the 
initial plane of the motion. 

After loosing orbital energy owing to the interaction with the en- 
vironment, the binary becomes hard when the semimajor axis 
reaches the value defined in Merritt et al. (2007), 


qd Fh 


KEZ 7 


an 
where q = Mj, y/Mj,,, = 0.03 is the SMBH mass ratio. In our 
model, r, ~ 156.6 pc, and therefore a, ~ 1.1 pc. Fig. 8 shows 
the distance between the two SMBHs and the semimajor axis of 
their orbit, together with the emitted GW power, as a function 
of time. In absolute values, the slopes of the three lines start to 
increase rapidly when the separation approximately reaches ap. 

Through the comparison with the case in which the calcula- 
tion does not take the PN terms into account, we can infer that 
most of the orbital decay is due to the dynamical friction. We 
determined that the only interaction with the background stars 
can carry the binary at a distance of a few times the character- 
istic Schwarzschild radius of the binary, which in our case is 
R, = 1.38 x 1075 pc. In the absence of PN terms, the binary 
would start to slowly shrink over a time of tens of million of 
years after this. The emission of gravitational waves rapidly re- 
duces the binary orbital energy and speeds up the decay: in a 
few thousand years, the distance plummets from few times R, 
to zero. According to our results, the merger between the two 
SMBHs occurs 16.6 Myr after the formation of Milkomeda, in 
the same range of timescales as was found by Khan et al. (2016) 
for similar SMBH mergers. 

The blue and red curves in Fig. 9 show the evolution of the 
eccentricity e and the inverse semi-major axis 1/a, respectively. 
The eccentricity starts at a value of 0.7 and drops to zero, as the 
binary shrinks and circularizes. 

As expected, when the binary becomes hard, the hardening rate, 
defined by Merritt et al. (2007) as 


dfi 
S= uNa) 


suddenly increases from ~ 4 pe™!Myr™! to ~ 4 x 104 pe~!Myr™! 
in the last phases before the merger. 

However, because the external environment in ARWV is not mod- 
eled with an ensemble of particles, we cannot reproduce the or- 
bital energy loss due to the stochastic gravitational encounters of 


(4) 
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Fig. 8. Evolution of the relative distance between the two SMBHs (blue 
line) and the semimajor axis a of their orbit (red line), together with 
the power emitted in the form of GWs (green line). The time is counted 
starting from the galaxy merger. 
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Fig. 9. Inverse of the semimajor axis of the orbit of the SMBHs (in pc~!) 
and its eccentricity as a function of time. The ordinate scale is the same 
for a and e. 


the binary with the field stars. Our estimation of s, and of |de/dt| 
as well, in the phases soon after the binary becomes hard is there- 
fore underestimated. We can obtain an estimate of the hardening 
rate due to the energy exchange during the encounters, in the as- 
sumption of a background with a fixed and uniform density p and 
uniform velocity v, through the approximated formula (Quinlan 
1996; Gualandris et al. 2016) 


G 
ss Ey 
v 


(5) 


where H is a dimensionless hardening coefficient with a nearly 
constant value of ~ 16 for hard binaries. In our case, this is 
s ~ 0.21 pc-!Myr“|, in accordance with those obtained in simi- 
lar simulations by Khan et al. (2018). 


5.2.2. Gravitational wave emission 


The power emitted in the form of GWs as a function of time, 
shown with the green line in Fig. 8, was obtained as the energy 
loss rate, averaged over one orbital period, according to the for- 
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mula (16) of Peters & Mathews (1963), 


7 7 
a2 e). 


P) = 
sP) * 54° * 96 


(6) 
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The emitted power progressively increases when the semi- 
major axis drops below a, and reaches a maximum of about 
10” Lo just before the merger. 
Through a numerical integration of the Peters formula, we ob- 
tained the amount of energy that is radiated away during the pro- 
cess. In Fig. 10 we compare the results of this integration with 
the energy loss calculated by the ARWV code through the PN ap- 
proximation. The two curves agree well, with a fractionary log- 
arithmic variation below ~ 13%. The amount of energy emitted 
in the last phases of process is on the order of 10* J. 

We repeated the simulation with the ARWV code for two ex- 
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Fig. 10. Energy radiated during the last phases of the SMBH merger. 
The green line shows the mean emitted energy, obtained from the inte- 
gration of Eq. (6), and red points show the emitted energy as output of 
the ARWV code. 


treme configurations of the SMBH spin vectors: one for paral- 
lel spins, and the other for antiparallel spins. As expected, the 
time of the merger is not affected by the spin orientation, but we 
found that the final merger remnant gains a different recoil ve- 
locity. In the case of parallel spins, the recoil velocity of the rem- 
nant is v, = 2.2 km/s, while for antiparallel spins, it is v, = 24.8 
km/s. The magnitude of the recoil velocity in both cases is small, 
mainly because of the low SMBH mass ratio (Healy & Lousto 
2018; Chassonnery & Capuzzo-Dolcetta 2020). The final SMBH 
is thus unable to escape from the center of the galaxy because 
the central escape velocity for our model is ve = 3.7 x 10° km/s. 
Therefore, the giant elliptical galaxy Milkomeda will continue 
to host an SMBH in its center, as obtained also by Arca Sedda 
et al. (2019). 

We also investigated the possibility of observing a GW signal 
from a merger between similar SMBHs in the near Universe. The 
frequency-characteristic strain evolution is shown in Fig. 11 and 
overlaps the sensitivity curve of different GW detectors, such as 
the Pulsar Timing Array (PTA, Hobbs et al. (2010)), the Square 
Kilometer Array (SKA, Johnston et al. (2007)), the Laser Inter- 
ferometer Space Antenna (LISA, Amaro-Seoane et al. (2017)), 
the Deci-Hertz Interferometer Gravitational wave Observatory 
(DECIGO, Kawamura et al. (2011)), and the wAres microhertz 
detector (Sesana et al. 2019). A comparison between our mod- 
eled signal and the detector sensitivities shows that mergers sim- 
ilar to the one we expect to witness in Milkomeda can be bright 


sources in ground-based detectors such as the PTA, or in the next 
decade, the SKA, provided that they take place roughly within 
1 Mpc. However, farther away in space, it becomes clear that 
the only possibility of observing this type of SMBH mergers is 
using space-borne detectors. Mergers that occur up to redshift 
z < 2 might be observable by LISA, if with an allegedly low 
signal-to-noise ratio, but they might shine bright to a micro-hertz 
observatory such as the Ares detector concept design. 
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Fig. 11. Characteristic strain of a GW signal emitted by a similar SMBH 
merging binary as a function of frequency for four different redshifts. 
The GW signal is computed starting from a semimajor axis of aọ = 1980 
AU and an eccentricity of eọ = 0.064, corresponding to a time of t = 4 
yr before the merger. The sensibility curves of the main GW detectors 
are also shown. 


6. Discussion and conclusions 


We studied the future evolution of the system composed of the 
Milky Way and M31 (Andromeda galaxy) in relation with the 
gravitational effects due to the intergalactic background. This 
can shed light not only on the forthcoming dynamics of our spe- 
cific galaxy system, but also in general on the correlation be- 
tween the galaxy interaction timing and the environmental prop- 
erties, the galactic structure, and the initial conditions. We used 
the high-precision N-body HiGPUs code, properly modified to 
account for the friction exerted on the bodies in the galaxies by 
the diffuse background, to calculate the evolution of the galac- 
tic hosts and their central SMBHs until the merger of the two 
galaxies. At this stage, we used the galaxy density, the veloc- 
ity distribution, and the SMBH orbital parameters for a further 
dynamical simulation. With the ARWV code in its most recent 
version, which considers PN treatment and the BH spins, we 
reproduced the orbital decay of the SMBH binary in the post- 
merger galactic background. The aim was not simply to estimate 
the likelihood and future time of an eventual merger, but also to 
determine the fate of the two massive black holes hosted at the 
two galactic centers, whose estimated mass is 4.31 x 10° Mo in 
our Galaxy and 1.4 x 108 Mo in M31, giving a mass ratio 0.03. 
We summarize our results below. 


1. The time evolution of the MW and M31 orbits is such that 
the first close approach of the two galaxies will occur in 4—5 
Gyr, with a weak dependence on the characteristics assumed 
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for the background density, the dimension of the halos, and 
the initial velocity. 

2. The time for the completion of the two galaxy merger in- 
creases significantly with the relative velocity transverse 
component, which is an ill-determined observational quan- 
tity. According to the most recent estimates, however we can 
conclude that the MW and M31 will merge in ~ 10 Gyr. 

3. The dimensions of the galactic halos play an important role 
in the merger time: larger halos cause a significant orbital 
energy dissipation and accelerate the decay, at least until the 
distance between the two galaxies is on the same order as the 
size of the halos. 

4. As expected, due to the collisionless nature of the encounter 
and merger, the overall post-merger density profile is not 
very different from a mere mass average of the two profiles 
of the MW and M31. 

5. After the merger, the two SMBHs were left orbiting on a 
mutual orbit of eccentricity e ~ 0.7 and semimajor axis a ~ 
160 pc, which stalled because the resolution of the N-body 
simulation is insufficient. 

6. The following fate of the SMBH pair was followed by a PN 
simulation that showed how efficiently (in less than 17 Myr) 
dynamical friction braking leads the two SMBHs to the so- 
called hard binary phase, when subsequent orbital decay is 
given by energy dissipation by GW emission down to the 
final merger and recoil kick. 

7. When we also considered antiparallel spins for the SMBHs, 
the recoil kick velocity was below 25 km/s (two orders of 
magnitude lower than the central escape velocity), which 
leaves the BH remnant confined in the inner potential well 
of the galaxy. 

8. Types of SMBH orbital decays similar to those studied here 
show a very high power of GW emission because of the high 
masses of the BH involved. This high power is shifted toward 
very low frequency. This GW emission would in principle 
be observable only with future GW ground-based detectors 
such as the PTA and the SKA or with space interferometers 
such as LISA, but the redshift range for the detection should 
bel <z<2. 


On the basis of our results, we are now able to determine 
the most feasible scenario of the future of our own Galaxy and 
its central SMBH. Our new estimate of the time required for the 
completion of the MW-M31 merger means that the life of the 
Local Group is slightly longer than previously believed. A final 
result is that unfortunately, the Sun will not live long enough to 
witness the formation of Milkomeda and will therefore not be 
part of the new galaxy. 
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